Quantifying uniaxial prestress and waveguide effects on dynamic elastography estimates for a cylindrical rod

Dynamic elastography attempts to reconstruct quantitative maps of the viscoelastic properties of materials by noninvasively measuring mechanical wave motion in them. The target motion is typically transversely-polarized relative to the wave propagation direction, such as bulk shear wave motion. In addition to neglecting waveguide effects caused by small lengths in one dimension or more, many reconstruction strategies also ignore nonzero, non-isotropic static preloads. Significant anisotropic prestress is inherent to the functional role of some biological materials of interest, which also are small in size relative to shear wavelengths in one or more dimensions. A cylindrically shaped polymer structure with isotropic material properties is statically elongated along its axis while its response to circumferentially-, axially-, and radially-polarized vibratory excitation is measured using optical or magnetic resonance elastography. Computational finite element simulations augment and aid in the interpretation of experimental measurements. We examine the interplay between uniaxial prestress and waveguide effects. A coordinate transformation approach previously used to simplify the reconstruction of un-prestressed transversely isotropic material properties based on elastography measurements is adapted with partial success to estimate material viscoelastic properties and prestress conditions without requiring advanced knowledge of either.


I. INTRODUCTION
Dynamic elastography methods-using noninvasive optical, ultrasonic, or magnetic resonance imaging modalities-aim to quantitatively map the shear viscoelastic properties of materials.These properties in biological tissues are often altered by disease and injury.When considering larger (relative to wavelength) regions of interest, like the liver or brain, boundary effects often can be ignored.However, as elastography expands to other anatomical regions where dimensions in at least one direction are smaller or of comparable length to bulk shear wavelengths-such as in slender skeletal muscles, blood vessels, the heart wall, and the cornea-boundary effects become non-negligible and must be considered.Researchers using optical elastography to assess the viscoelastic properties of the cornea have long recognized this, adapting models to include waveguides by treating the cornea as a plate-like structure.3][4] Limited studies on cardiac elastography have also acknowledged the frequency-dependent (i.e., wavelengthdependent) waveguide behavior of the heart wall. 5ften, when elastography studies are done under varying nonzero quasi-static prestress conditions, observed changes in mechanical wave behavior are attributed solely to the nonlinear property of the tissue: it is observed that its shear and viscous constants are highly dependent on the static load and associated deformation.A recent article provides a summary of the literature relevant to this issue, in particular for uniaxially prestressed cylindrically-shaped structures, as well as biaxially prestressed plate-like structures. 6In another recent article, the impact of neglecting the prestress effect in cornea elastography is quantified. 79][20] We investigate the confounding effects of finite dimensions and prestress on elastography measurements.We articulate and evaluate a strategy for decoupling prestress and waveguide effects from estimates of material shear viscoelastic properties based on circumferentially-, axially-, and transverselypolarized wave motion in the cylinder.

II. THEORY A. Mechanical wave motion in a uniaxially prestressed linear viscoelastic material
Most dynamic elastography methods assume that the measured transverse wave speed or wavelength for small amplitude (linear theory assumption) motion is directly a) Email: troyston@uic.edurelated to the material's elastic or viscoelastic properties. 21ssuming isotropy, homogeneity, and neglecting boundary effects or variation in density q in a nearly incompressible viscoelastic material, the frequency-dependent shear wave phase speed c x ½ for harmonic excitation at circular frequency, x, is Here, k sh x ½ is the complex-valued, frequency-dependent shear wave number, and l x ½ is the complex-valued, frequency-dependent shear modulus, comprised of the shear storage modulus, l R x ½ , and the shear loss modulus, l I x ½ , such that l x . The attenuation rate of the wave as it propagates is governed by the imaginary part of the wavenumber: Imag k sh ½ .In a viscoelastic material, both the shear storage and loss moduli affect both the phase speed and attenuation rate.In a purely elastic material, l I ¼ 0, there is no attenuation and the phase speed is independent of frequency (nondispersive) and reduces to c ¼ ffiffiffiffiffiffiffiffi l=q p .While some linear studies have assumed pure elasticity (no viscosity), often their analyses are generalizable to the linear viscoelastic problem for harmonic motion by simply adding the imaginary shear loss modulus to form the complex shear modulus.This approach is used in the present analysis.
Consider the introduction of a uniaxial static prestress r parallel to the z-axis (Fig. 1).If the static deformation due to the prestress is assumed sufficiently small such that higher order nonlinear terms can be neglected, the governing equations expressed in polar coordinates r; u; z ð Þ are as follows: where u; v; w ð Þ denote incremental displacements in the r; u; z ð Þ directions, S rr ; S uu and S zz are normal stresses in the respective directions, and S ru ; S uz ; S zr ; S rz ; S zu and S ur are shear stresses with the first subscripted dimension referring to the shear surface and the second subscripted dimension denoting the shear direction.Finally, subscripted r; u; z; t ð Þ after a comma refer to partial derivatives with respect to that spatial or time dimension, 22 S rr;r þ 1 r S rr À S uu þ S ur;u ð Þ þ S zr;z ¼ qu ;tt ; (2) The normal and shear stress terms are as follows: Here, k is the volume elasticity of the material, a Lame constant.It is closely related to the material's bulk modulus j ¼ k þ ð2=3Þl.In biological soft tissue or other "nearly incompressible" materials, k and j are multiple orders of magnitude greater than l such that the Poisson's ratio for the material approaches, but does not equal 0.5.The terms l r and k r are coefficients that regulate material behavior in the presence of an initial stress, which in this case is uniaxial prestress r aligned with the z (cylinder) axis, as shown in Fig. 1.Specifically, they account for a linear dependence of l and k on r.In the Refs.22 and 23, these terms are b 1 and b 2 , and b 1 is related to A used in other references, [24][25][26] where its value in soft tissue-like nearly incompressible materials can range between negative and positive values.It has been hypothesized that its value depends on microstructure and thus can reveal material changes not captured by l. 24 As noted, the previous formulation neglects higher (3rd, 4th, etc.) order nonlinear terms in the strain energy function and the limitations of this simplification will be evident in the numerical and experimental studies of Secs.III and IV.Incorporation of the higher order terms, which have been detailed in previous studies investigating bulk wave motion without waveguide effects, 23 is left for future study.][26] In the present study, we will assume harmonic steady state motion.By using the separation of variables one form of the solution for u; v; and w, the symmetric form, is u r; u; z; t The anti-symmetric form is u r; u; z; t In the following subsections, we consider specific cases of torsionally-, axially-, and transversely-polarized harmonic wave excitation that leads to simplifications of the previous expressions.

B. Torsionally-polarized wave front
Following Section 8.2.2 of Graff, 27 consider waves that are torsionally-polarized and axisymmetric, meaning no variation with u.The anti-symmetric Eqs. ( 17)- (19) reduce to the case of n ¼ 0, and thus u r; u; z; t where J 1 denotes a Bessel function of the first kind of order 1.The radius of the cylinder will determine allowable values of b > 0, for which there will be an infinite number such that, in general, v r; z; t This is satisfied for the case that b ¼ 0 and also imposes that bRJ 0 bR ½ ¼ J 1 bR ½ , which leads to values b 1 R ¼ 5:136; b 2 R ¼ 8:417; b 3 R ¼ 11:62; … While there are an infinite number of solution sets, given attenuation due to viscosity we expect the lowest wavenumber values, associated with the longest wavelengths, to dominate as we move further from the source of excitation.Additionally, the geometric configuration of the source of the excitation may preferentially drive certain wavenumber solutions.
With this form for v r; z; t ½ , Eq. ( 3) reduces to the following: Here, n and b denote complex wavenumbers, with n in the z-direction, and b in the r-direction.Equation (22) shows how they are related to the shear wave number k sh .For the case that b ¼ 0 and r ¼ 0 we have n 0 ¼ k sh .In other words, the torsional waves propagate along the cylinder axis exactly as bulk shear waves, with wavelength k sh equal to 2p=Real k sh ½ .If there is an axial prestress r then the wavelength will equal Assuming the material to be nearly incompressible, axial prestress r and axial prestrain are related by r ¼ 3l 0 , where l 0 is the static (real number) value of l x ½ at x ¼ 0. If one is able to acquire estimates of n 0 based on fitting e Àjn 0 z to the measured wave profile along the cylinder axis z at a given r and u at both unstressed and prestressed conditions at multiple frequencies, then l x ½ and r can be independently determined, even without being able to directly measure r under static load conditions.(A direct measurement of r may not be possible in vivo or in any situation where the cylindrical structure is connected at its ends to another structure and it is impossible to measure the stress condition between them.)Specifically, by selectively driving (or measuring via appropriate filtering) only the b ¼ 0 mode, measurements at multiple frequencies without prestress enable a determination of l x ½ , including (possibly by extrapolation with an assumed rheological model) l 0 ¼ l x ¼ 0 ½ .Then, we can use the fact that r ¼ 3l 0 and Eq.(23) to independently determine values for r and l r , assuming we can estimate from imaging.Numerical and experimental studies in the Secs.III and IV, respectively, illustrate and evaluate this process.

C. Axially-polarized axisymmetric wave front
We next consider the case of an axisymmetric axiallypolarized wave front, which theoretically could be driven by axial oscillation of the rigid cuff encircling the cylinder over a finite axial length.This again results in n ¼ 0, but now with the symmetric form, Eqs. ( 14)- (16), which reduce to u r; z; t w r; z; t Following Graff, 27 Section 8.2.2, the boundary condition at the outer radius R means that S rz ¼ S rz ¼ 0, which is satisfied by the following forms for U r ½ and W r ½ : Here again, there will be certain allowable values of a 1 ; a 2 ; …; a 1 , b 1 ; b 2 ; …; b 1 , and n 1 ; n 2 ; …; n 1 .Inserting these forms for u and w into Eq.( 4) results in the following, where Bessel functions with a and b are separated into two separate equations, Eqs. ( 28) and (29), that both must be satisfied, Equation (28) shows the coupled relationship between axial wavenumber n, radial wavenumber b, and shear wavenumber k sh .Equation (29) relates n, another axial wavenumber a; and compression (longitudinal) wavenumber k p .If prestress r ¼ 0, the previous equations reduce to expressions found in Graff. 27nlike in the torsional case, even with selective excitation of, say, the lowest values for a 1 , b 1 ; and n 1 , there is not a case with a or b ¼ 0. Nonetheless, Eqs.(28) and (29) do show the relationship between shear, compression, axial, and radial wavenumbers, and how this is effected by uniaxial prestress.Focusing on Eq. ( 28), if values for b 1 and n 1 can be determined at multiple frequencies with no prestress present, they can be used to then estimate l x ½ .Then, prestress is added to determine r ¼ 3l 0 based on measuring .Following that, Eq. ( 28) can be used at different prestress levels to estimate l r and k r .Numerical and experimental studies in the following Secs.III and IV, respectively, illustrate and evaluate this process.

D. Transversely-polarized axisymmetric wave front
For wave symmetric to the x-axis in Fig. 1, we restrict ourselves to the case of n ¼ 1, and Eqs. ( 14)-( 16) reduce to the following: Unfortunately, this case does not simplify to the same degree as did the torsionally-and axially-polarized cases.All three displacements u, v, and w, are present, and the resulting forms for U r ½ , V r ½ , and W r ½ remain more complex. 27Further analytical development of these solution forms in the presence of uniaxial prestress is left for future study.In the present article, we instead consider simplifying theories that may yield insight as the cylinder's radius becomes small as compared to the wavelengths being considered.We specifically review "thin" and "thick" beam theories under uniaxial prestress conditions.
1. Prestressed one-dimensional "thin" waveguide under transverse excitation By assuming that the radius R of the cylinder is small enough such that bR ( 1, then so-called one-dimensional beam theory approximations can be applied.The simplest of these is the Euler-Bernoulli thin beam theory, which allows for the incorporation of prestress r.Referring again to Fig. 1, the pre-tensioned Euler-Bernoulli thin beam theory described in Section 3.3.4 of Graff 27 for x-polarized transverse wave propagation of the beam along its z-axis is the following: Here, I is the area moment of inertia about the y-axis ½I ¼ ðp=4ÞR 4 for a circular cross section of radius R], and A is the cross-sectional area in the x-y plane.The general solution form is: v ¼ Ve jðxtÀnzÞ where n has four possible solutions, We have two propagating waves (a) in the þ or -z-direction, and two non-propagating (near field or evanescent) waves (b) in the þ or -z direction.For the propagating waves, the phase speed will be: c ph ¼ ðx=Real a ½ Þ. Taking the limit that r=l ( 3R 2 =32 we see that a ¼ x=a ð Þ 1=2 ; and thus for the elastic case, c ph ¼ ðx 1=2 =2Þð3lR 2 =4qÞ 1=4 , which is the classic thin (Euler-Bernoulli) beam transverse vibration solution.On the other hand, taking the limit of tension r=l ) 3R 2 =32 we then drop EIu ;zzzz from Eq. ( 33) and reformulate the solution to find that there are two propagating solutions with c ph ¼ r=q ð Þ 1=2 .This is the classic transverse thin string vibration solution.A case between the extremes of either neglecting rA or EI still does not match the phase speed of bulk shear waves, given by Eq. (1).

Prestressed one-dimensional "thick" waveguide under transverse excitation
The thin beam formulation presented in Sec.II D 1 is only a reasonable approximation when the wavelength of the propagating transverse waves is at least an order of magnitude greater than the beam's cross-sectional radius R, or equivalent radius for a noncircular cross section.A formulation yielding a valid approximation for shorter transverse waves is based on the Timoshenko beam theory, which allows for shear deformation and accounts for rotational inertia.Incorporating prestress into the formulation for Timoshenko beam theory given in Section 3.4 of Graff, 27 we have the following: Here, v is the Timoshenko shear coefficient, equal to 10=9 for a circular cross section.Consider harmonic motion at frequency x and the general solution form is u ¼ Ue jðxtÀnxÞ , and n has four possible solutions, In this case, we cannot separate the wavenumber solutions into propagating and non-propagating pairs.Rather, there will be certain prestress-dependent frequency ranges where we have only one propagating pair and others where we have two propagating pairs.
E. Accounting for prestress using "transformation acousto-elastography" (TAE) Equations ( 23), (28), and ( 29) "hint" at an alternative strategy to solve the inverse problem of independently identifying inherent material shear viscoelastic properties and prestress based on measurements of wave motion.In previous studies on transversely isotropic materials not under prestress, [28][29][30][31][32] the last author has shown that, by distorting the geometry based on direction and polarization-dependent planar phase speeds one can then solve an equivalent isotropic problem.This approach to the anisotropic problem is "Transformation Elastography."Uniaxial prestress causes a similar, though not identical, direction and polarization dependence of the planar shear wave phase speed.The same approach is adapted to the acoustoelastic problem here, which also has the added complexity of waveguide effects, and is evaluated in numerical and experimental studies noted in the following.
The analysis in the previous sections shows that the case of torsionally-polarized waves may be the simplest and is a good place to start.Considering that the material is nearly incompressible, the prestress results in a static strain of the cylindrical phantom of unstressed length L and radius R, changing its axial length to L 1 þ ð Þ¼ Lð1 þ ðr=3l 0 ÞÞ and its radius to By driving or measuring only b ¼ 0 torsional wave motion, if one then distorts the axial length by dividing it by Real ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi , then fits a solution to e Àjn Ã 0 z in order to identify n Ã 0 , where an asterisk * is added to denote this is based on a distorted geometry, then we see that n Ã 0 ¼ k sh , where k sh was obtained under the unstressed condition.From the degree of distortion needed to achieve this equality, we directly determine l r .
Note, l x ½ (l 0 ; l a ; a) needs to be determined for the no prestrain ( ¼ 0) case for multiple frequencies.So, instead of adjusting the geometry distortion by varying l r as we propose here to match with the unstressed case, another approach is to adjust l r to match the measured n 0 in the undistorted geometry using Eq. ( 23).This will lead to the same result as arrived at using the geometry distortion (TAE approach).However, we propose that the TAE concept, applicable to more than the torsional polarization case, may help conceptually as an analysis tool.
Next, consider the axisymmetric, axially-polarized wavefront.In this case, the geometry is distorted by dividing the axial length by Real½ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Assume that, within the region bounded by the rigid cuff, we have purely radially converging axially-polarized waves that can be fit to J 0 b Ã r ½ in order to solve for b Ã .This follows from Eq. ( 27) since it is assumed that n Ã ¼ 0. Additionally, with this assumption, we see from Eq. ( 37) that b Ã ¼ k sh : Thus, like in the torsionally-polarized case, we can match b Ã to the value of k sh that was obtained with no prestress by adjusting l r .This becomes another means of identifying l r .Then, by measurement of the wave field generated in both the radial and axial directions with and without axial prestress present, we can use Eqs.( 26)-( 29) to estimate k r .However, there are an infinite number of solutions to this transcendental set of equations.Numerical and experimental studies in Sections III and IV evaluate whether this approach can be used to determine k r .Furthermore, while a similar analysis has not been completed for the case of transverselypolarized (flexural) waves, numerical studies detailed in the following provide insight into the relationship between prestress and l r and k r in this case.

A. Methods
An analytical and numerical case study was conducted to understand the interactions between uniaxial prestress and waveguide behavior, as well as to evaluate the TAE approach introduced in Sec.II E. Wave motion in the cylindrical phantom shown in Fig. 1, with parameter values provided in Table I, was simulated with a numerical finite element (FE) approach using ANSYS Mechanical APDL Version 2022 R1 (Ansys, Canonsburg, PA).These parameter values were chosen to be similar (though not exactly matching) to those in experimental studies described in Sec.IV.For simulation of torsionally-and axially-polarized harmonic excitation, an axisymmetric mixed u-P formulation was used with Plane183 8-node elements with individual element side dimensions of 0.5 mm.(Both displacements and hydrostatic pressure are taken as primary unknowns in the mixed u-P formulation, which is recommended for nearly incompressible materials.)For transverse x-polarized excitation [referring to Fig. 1(c)] a mixed u-P formulation was used with Solid273 8-node by six circumferential plane generalized axisymmetric elements with individual element side dimensions in the axial and radial direction of 0.5 mm.In all cases, first, the response to the static preload is solved by accounting for geometric nonlinearity (nlgeom, on).One end of the cylinder was fixed in the axial (z)direction and the other end was incrementally displaced in the z-direction, solving the problem in steps until the desired end displacement was reached that resulted in a uniform uniaxial strain throughout the model.This was done using the Gent model 33 defined by, j; l 0 ; and J m , or using linear elastic properties, E; ; and l 0 , where Young's modulus E and Poisson's ratio were chosen to be consistent with the value of j used for the Gent model.For axial strains up to 20%, differences between the Gent and linear model, as quantified by different predicted axial stress values, were about half the percentage of the strain.Specifically, the axial stress predicted by the Gent model exceeded the stress predicted by the linear model by 1.25% for 2.5% axial prestrain, 2.54% for 5% prestrain, and 11.2% for 20% prestrain.While the Gent model is more accurate, it was necessary to use the linear model in Ansys for the type of linear harmonic perturbation analysis conducted here, which required that the shear modulus value be changed between the static and harmonic analysis phases.This was needed due to the frequency-dependence of the shear modulus.
Specifically, once the static analysis was done, the solution routine was exited and then re-entered, recovering the prestrained geometry and a "total" stiffness matrix that is the sum of the initial stiffness matrix due to material properties plus terms added due to the prestress condition. 34Then, the shear and Young's modulus values are updated, which in turn updates the "total" stiffness matrix.This is done before turning geometric nonlinearity off (nlgeom,off) and entering the harmonic solution routine using this modified total stiffness matrix acquired at the end of the static solution routine.Now, harmonic u; z; or x-direction displacements were applied.For torsional and transverse (flexural) excitation, one end of the phantom was specified to have upolarized or x-polarixed motion, respectively, in order to preferentially drive the b ¼ 0 torsional mode or n ¼ 1 flexural mode, respectively.For axially polarized excitation, nodes on the outer surface of the phantom representing a rigid axially-oscillating cuff in contact with the phantom from 35 to 50 mm axially were given harmonic z-polarized displacement inputs.In addition to these simulations, we also computed the response to an axially-and to a transversely-polarized 10 mm long line segment displacement located at the geometric center of the cylinder, as an approximation of dynamic elastography using modulated radiation force of ultrasound to create a force vector at the focal region of the ultrasound probe.
For the harmonic analysis, a "Fractional Voigt" rheological model is assumed 35 and the shear modulus is modified from l 0 to l x ½ ¼ l 0 þ l a x a cos pa=2 ½ and viscous (beta) damping is specified by taking the ratio of the imaginary to the real part of the shear modulus and dividing by x, the harmonic circular driving frequency.Specifically, the beta damping value is B. Results The FEA-calculated in-phase steady-state torsionally (u)-, axially (z)-, and transversely (x)-polarized response over a sagittal slice (FOVs in Fig. 1) is shown in Figs. 2, 3,  and 4, respectively, at 0%, 10%, and 20% axial strain levels at 600 Hz for the 35 mm diameter phantom.
The acoustoelastography challenge is to estimate the complex frequency-dependent shear modulus and the prestress and prestress-dependent coefficients, l r and k r , based only on assuming a known density q and using the complex wave images at multiple frequencies and prestress levels, from which we can also determine axial prestrain .The procedures described in Sec.II E are used to estimate the rheological properties first at multiple frequencies under zero prestrain in order to obtain estimates of l 0 ; l a ; and a, and then with prestrain, in order to estimate prestress r, and coefficients l r and k r .We start with a curve fit estimate of the wave profile along the z-direction during torsionallypolarized harmonic excitation, followed by fits of the wave profile along r-and z-directions during axially-polarized harmonic excitation.Transversely-polarized (x-direction) and line segment force studies following this are used to identify the effect of prestress r on flexural wave motion, and to assess the combined prestress and waveguide effect on shear wave motion generated from a focused source, as an approximation of a modulated radiation force of ultrasound source.Results at 600 Hz are summarized in Tables II and III for the 35 and 20 mm diameter phantoms, respectively.
Curve fits were done in MATLAB using the "lsqcurvefit" command.Fits for the torsionally-polarized excitation were along 5 cm axial lines, averaging fits from lines spaced 0.5 mm apart from r ¼0.5 mm to r ¼ R. First, with no prestrain, it was confirmed that this fitting approach conducted at multiple frequencies enabled the identification of l 0 ; l a ; and a of the Fractional Voigt rheological model.In Fig. 5, the black dashed line shows the real (storage) part l R versus the imaginary (loss) part l I of the shear modulus as frequency is varied.The blue x's denote the estimated values of l R verus l I based on curvefitting FE simulations at 500, 600, and 700 Hz for the 35 mm diameter phantom (they lie on the black dashed line as expected).By drawing a line through these x's, the intercept with the horizontal axis is l 0 .The real and imaginary values at a given frequency can then be used to determine l a and a.
For the nonzero prestrain cases, the length distortion was iterated by adjusting l r as discussed in Sec.II E until n Ã 0 at 2.5% prestrain matched k sh (previously obtained with 0% prestrain) minimizing the least square error.It was found that setting, l r ¼ À0:75, gave the best fit, with a difference in predicted k sh of less than 0.1% as compared to the unstressed case, for both the 35 and 20 mm diameter phantoms, as provided in Tables II and III, respectively.This close agreement was also present for 5% prestrain, but results diverged to about 0.4% and 1.4% difference at 10% and 20% prestrain, respectively, possibly due to geometric nonlinearities that are not accounted for in the linear-based approach.Figure 5 also shows how corresponding estimates of l change as prestrain is increased to 5% without the TAE correction (red circle) and with the TAE correction (green asterisks).
For the axially-polarized case, first a line fit was done along the radial line within the excitation zone (cuff).It was found that l r ¼ À0:75 also gave a good fit, with a difference in predicted k sh of less than 0.4% as compared to the unstressed case, for both the 20 and 35 mm diameter phantoms tested at 2.5% and 5% prestrain.This provided confidence in the theoretical developments of Sec.II.Results diverged up to 1.2% difference at 10% and 20% prestrain.
Next, averaged axial line fits were located at 0.5 mm spacing from r ¼ 0 to r ¼ R extending 30 to 50 mm distance from the excitation cuff (based on the unstrained geometry) in an attempt to reduce near field effects closer to the source, which may contain larger contributions from n > 1 of the a 1 ; a 2 ; …; a 1 , b 1 ; b 2 ; …; b 1 , and n 1 ; n 2 ; …; n 1 solutions.Then, we averaged fits along radial lines 0 < r < R taken at axial distances from 30 to 50 mm from the source in 0.5 mm increments.With l r already determined, the length distortion was iterated solely by adjusting k r until at 2.5% prestrain the values of b Ã and n Ã satisfied Eq. ( 37) where the value used for k sh was based on 0% prestrain.As expected, for this case estimates of k sh under unstressed conditions r ¼ 0 ð Þ were not accurate, being 5% too high and 7% too low for the 35 and 20 mm diameter phantoms, respectively.However, by setting k r ¼ 2:0 it was found that cases with 2.5% prestrain yielded predictions of k sh consistent with the unstressed case within 0.2%.This difference increased to about 6% at 20% prestrain.
For the transversely-polarized (flexural wave) case, fits along axial lines were located in 0.5 mm increments from r ¼ 0 to r ¼ R extending 30 to 50 mm distance from the excitation (based on the unstrained geometry), as in the axially-polarized study, in order to estimate n 1 .The same length distortion used in the axial study, Real½ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi , was applied here in order to estimate n Ã 1 .As in the axially polarized case it is not surprising that the values of n 1 obtained under zero prestress do not align with k sh .They are influenced by waveguide effects (small diameter of the phantom relative to shear wavelength).However, as in the axially-polarized case, by setting k r ¼ 2:0 it was found that cases with 2.5% prestrain yielded predictions of k sh consistent with the unstressed case within 0.4%.This difference increased to only about 1.5% at 20% prestrain.This suggests that the transverse waves are affected by uniaxial prestress in the same way as the axially-polarized waves, and thus may provide a more direct estimate of k r .
Theoretical estimates of n using "thick" beam theory were within 10% of values found using the finite element analysis (FEA) and reduced in value with respect to increasing prestress, like in the FEA simulations, though to not the same degree observed in the simulations.As expected, thick beam theory was a better match for the 20 mm versus the 35 mm diameter phantom.Thin beam theory results were poor in both cases.This is expected since the wavelengths generated at 600 Hz were comparable to radial dimensions.They would need to be an order of magnitude greater than radial dimensions for thin beam theory to be a good approximation.
To approximate focused modulated radiation force of ultrasound excitation, harmonic 10 mm long line displacements were input at the geometric center of the cylinder in axially-and transversely-polarized directions, as noted in Fig. 1.The response along a line orthogonal to the truncated line source was then fit to determine the complex wavenumber.If we assume that the truncated line source is an infinite  plane in an infinite medium (neglecting finite boundaries) and is creating a planar shear wave front the axiallypolarized wave traveling in the radial direction with wavenumber b should be governed by In other words, b . This is in agreement with Eq. ( 28) for the case that n ¼ 0. Values reported in Tables II and III show that, indeed, b Ã at 2.5% prestrain is consistent with the value of b determined at 0% prestress, but these values differ from k sh by about 5% in terms of the real part, which is inversely proportional to shear wavelength 2p=Real k sh ½ ð Þ .The imaginary value for b is an order of magnitude greater than k sh .This is expected since, not only is there attenuation due to viscosity, but also due to geometric dispersion as the finite source spreads out in three dimensions, including the axial direction, as it propagates.
Likewise, if we assume that the truncated line source is an infinite plane in an infinite medium (neglecting finite boundaries) and is creating a planar shear wave front the radially-polarized wave traveling in the axial direction (direction of uniaxial prestress) with wavenumber n should be governed by In other words . This is not in agreement with Eq. ( 28) for the case that b ¼ 0, where in that case . Both n Ã a and n Ã b are provided in Tables II and III.We expect the imaginary parts of the wave numbers to be greater than that of k sh , which is 54.9 at 600 Hz, since there will be geometric spreading.However, in this case, the geometric spreading is more confined by the small diameter of the phantom relative to its axial length.Imaginary values are larger than k sh , but by not nearly as much as they were for the other line segment case.And, as expected since there is less spreading, the imaginary part of the wavenumber is less for the 20 mm versus the 35 mm diameter phantom, since there is less geometric dispersion ($120 versus $180 m À1 ).By looking at the trend of the real part of n Ã a and n Ã b as prestress increases, it appears that n Ã b is the better estimate for k sh , suggesting that finite diameter (waveguide effects) have impacted the transverse wave motion.

IV. EXPERIMENT A. Optical elastography methods
The experimental configuration, shown in Fig. 6, enables excitation and measurement of torsionally-polarized wave motion on the surface of the cylindrical phantom that is simultaneously prestrained along its axis.The fixture parts were designed in Solidworks (Solidworks 2021) and threedimensionally (3D) printed using a fused filament extrusion printer (Prusa Mk3, PRUSA REF, Prague, Czech Republic) on polyethylene terephthalate glycol (PETG).The level of prestrain, introduced by hanging weight (container of water) from the lower end, was measured using a caliper.
A scanning laser Doppler vibrometer (SLDV; PSV-400, Polytec, Waldbronn, Germany), is used to measure motion on the phantom that is in the direction of the laser beam, as described in previous studies. 36,37The phantoms used in this study were made of Ecoflex TM 00-30 (Smooth-On, Inc., Macungie, PA) platinum-catalyzed silicone.Ecoflex silicone bases were combined in a 1 A:1B ratio and cured at room temperature in 3D-printed molds in a vacuum chamber (5305-1212, Thermo Scientific-Nalgene, Rochester, NY) to remove any bubbles and improve phantom homogeneity.
Our group has characterized the dynamic shear viscoelastic properties of Ecoflex materials under infinitesimal (0% pre-strain) conditions over a wide frequency range. 35,38revious publications focused on Ecoflex-10 TM , a similar polymer, with greater viscosity, as compared to Ecoflex-30.For harmonic excitation of both materials using small perturbations about the unstressed state, we've found that a "fractional Voigt" rheological model best describes the frequency-dependent shear storage and loss moduli of the material.The 3-parameter fractional Voigt model is comprised of a purely elastic element of strength l 0 in parallel with a fractional order springpot that is defined by a and l a such that, at frequency, x ¼ 2pf rad=s ð Þ, the shear storage l R and loss l I moduli are given by the following equations: From Ref. 35, it was found that for Ecoflex-10, the following values describe shear viscoelastic properties over the range from 200 Hz to 7.75 kHz: l 0 ¼ 13:3 kPa, l a ¼ 2 kPa Á s a , and a ¼ 1=3.For an Ecoflex-30 sample, measurements conducted in the same way as in the Ref. 24 over the range from 200 Hz to 1 kHz yielded: l 0 ¼ 27 kPa, l a ¼ 1:5 kPa Á s a , and a ¼ 0:3.Thus, while Ecoflex-30 is "stiffer" (higher l 0 ) under static conditions, due to its lower viscosity the magnitude of its complex shear modulus increases at a slower rate with frequency (a ¼ 0:3), as compared to that of Ecoflex-10 (a ¼ 1=3).Note, while these values for Ecoflex-30 were used in the numerical finite element study of Sec.III, the torsional experimental study described in this section while confirming the appropriateness of the Fractional Voigt model, yielded different parameter values.This could be due to variations in batches of Ecoflex-30, as well as the different experimental configuration.SLDV measurements were made over a grid in a ROI as indicated in Fig. 6.By driving the two piezoceramic stack actuators (P842.10,PI USA, Auburn, MA) in phase only torsional motion should be excited in the ideal case that the experimental setup is perfectly symmetric in all aspects, which is never the case.It is expected that some flexural (transverse) wave motion will also be excited.If measurements at equal radial distances from the central axis at the same axial (z) position, are subtracted from one another (difference) this should double the b ¼ 0 torsional mode and eliminate the n ¼ 1 flexural mode.If the measurements are added, it should eliminate the torsional mode and double the n ¼ 1 flexural mode.This was done before taking line profiles axially to determine n 0 (m À1 ) in the torsional case and n 1 (m À1 ) in the flexural case from the same measurement.

B. Optical elastography results
Optical elastography measurements, processed to amplify torsional motion, as described in Sec.IV A, are shown in Fig. 7 at 0%, 5%, and 10% axial prestrain for the 35 mm diameter phantom.(Image is flipped vertically relative to diagram and photo in Fig. 5.) Corresponding estimates of n 0 (m À1 ) in the torsional case are provided in Tables II and III for the 35 and 20 mm diameter phantoms under different prestrain conditions.Estimates of n 1 (m À1 ) in the flexural case were found to be unreliable indicating that the torsional motion was significantly greater in amplitude, as was intended.The case with no prestrain yielded an estimate of the shear wavenumber of k sh ¼ n Ã 0 ¼ 583 À j54:4 (m À1 ) for the 35 mm phantom, and k sh ¼ n Ã 0 ¼ 568 À j110 (m À1 ) for the 20 mm phantom, suggesting the 35 mm phantom material matched values used in the FEA simulation within 1%, but the 20 mm phantom was slightly stiffer by 3% and significantly more viscous than the material parameter values used in the computational (FEA) simulations, where k sh ¼ n Ã 0 ¼ 588 À j54:9.A more comprehensive study of the 20 mm phantom was undertaken by conducting the analysis over the frequency range from 200 to 600 Hz under 0%, 2.5%, 5%, 10%, and 20% prestrain levels in order to fully evaluate the material parameter identification study, starting with identifying the unstressed fractional Voigt model parameters.Figure 8 shows the estimated complex shear modulus l x ½ under 0% pre-strain (blue lines) based on the experimental estimate of n 0 x ½ ¼ k sh x ½ and l x ½ ¼ q x=k sh x ½ À Á 2 .The start to break down at frequencies above 500 Hz, resulting in large fluctuations in the estimated complex shear modulus.
Figure 8 also shows the experimental estimate of l x ½ when a 10% prestrain is imposed.As expected, the value of l x ½ is increased by the prestrain.Distorting the axial length as in the numerical study by dividing it by where l r ¼ À0:5 results in l Ã x ½ , denoted by the green line in Fig. 8, which more closely agrees with the value for l x ½ obtained under the no prestrain case.Table IV summarizes experimentally fit wavenumber estimates n 0 and TAE-adjusted estimates of k sh ¼ n Ã 0 , all based on l r ¼ À0:5 for the same prestrain levels used in the numerical study.
For the experimental studies, the value of l r ¼ À0:5 provided a better fit than the value of l r ¼ À0:75 used in the FEA studies.The value of l r ¼ À0:5 is consistent with the initially-developed theory of acoustoelasticity articulated by Biot in his seminal text. 39More recent studies have shown that Biot's theory is a special case of a broader theory that also incorporates higher order effects. 23Experimental results in Tables II and III show that higher order terms are necessary for accuracy as prestrain increases with the estimate of k sh ¼ n Ã 0 deviating by more than 6% at 20% prestrain.

C. Magnetic resonance elastography methods
Axially-and transversely-polarized wave propagation in the same uniaxially prestressed cylindrical isotropic phantoms used in the optical measurements was studied experimentally using magnetic resonance elastography.All experiments were conducted in a 30 cm horizontal bore Agilent 9.4 Tesla preclinical magnetic resonance imaging (MRI) system using a 120 mm I.D. 600 mTesla/meter maximum strength gradient set and a 72 mm I.D. linear radio frequency coil.
The experimental fixture to hold the phantom in the magnet and simultaneously apply a known static tensile preload while also delivering vibratory excitation via a rigid cuff around the phantom at its mid-length is shown in Fig. 9.The fixture parts were designed in Solidworks (Solidworks 2021) and 3D printed using a fused filament extrusion printer (Prusa Mk3, PRUSA REF) on PETG to limit interference with magnetic resonance (MR) signals and keep parts light enough to be moved efficiently by the piezo stack actuator that delivered the vibratory excitation.It is an advancement of a prior design. 40Tensile pre-loading of the cylinder is shown in green in Fig. 9.Samples are gripped by two clamp types: a fixed clamp at the proximal end, and the tensioner clamp at the distal end.Both clamp types are slotted to allow wooden skewers to penetrate the sample, providing reasonably distributed tension, rather than surface tension that would be only available with typical clamps.The distal clamp is then attached to a reinforced nylon wire that is fed out the back of the magnet bore and attached to a pulley system where adjustable weight applies appropriate tension to achieve a desired measurable strain.
Simultaneous harmonic actuation (orange in Fig. 9) is achieved through a piezo stack (P842.10,PI USA LP) located proximally to the entry of the bore to limit interference from wiring, and attached to a Delrin countermass to maximize transference of motion to the phantom.The phantom is placed in a cylindrical tube and the actuator is attached directly to the piezo and slotted for adjustable cuffs (yellow) that uniformly and circumferentially grip the  sample, providing a fixed location of the wavesource that contacts the phantom over 15 mm in length axially.Harmonic vibration is produced along the same axis as the pre-tension.Note, however, that because of gravity and imperfections in this setup flexural (transverse) motion will also be excited.SLIM MRE 41 was used to acquire vibratory motion encoded in three orthogonal directions simultaneously.Sequence parameters were as follows: TR/TE 1600/16 ms; eight time steps were captured at even intervals over the period of the vibration frequency.A 250 mT/m motion encoding gradient (MEG), which was matched to the vibration frequency, was used with ten cycles.Measurements were repeated with inverted polarity gradients to subtract static field inhomogeneities.The data matrix size of 64 Â 64 Â 40 with an isotropic voxel size of 0.75 mm resulted in a field of view (FOV) of 48 mm Â 48 mm Â 30 mm.The piezo stack provided an input axial harmonic motion of $10 lm peak amplitude at 600 Hz.

D. MRE results
Sagittal views obtained using MRE of the axiallypolarized vibratory motion are provided in Fig. 10 for 0%, 5%, and 10% axial prestrain for the 35 mm diameter phantom.The field of view is the same as in the numerical study of Sec.III (except shortened to 3 cm along the axis).The images were analyzed using the same TAE approach detailed in Sec.III for the numerical study, with calculation results provided in Tables II and III.Estimates of b based on the radially converging axially-polarized wavefront within the cuff yielded results similar in value to the torsional studies for the 35 mm phantom, given that k sh ¼ b Ã .Like in the numerical studies, as prestrain is increased we find that the measured b value increases, as expected per equation (37).A corrected estimate of b Ã can be obtained by dividing the Corrected" values b Ã for nonzero prestrain cases more closely match the zero prestrain case, as can be seen in the last row of Tables II and III.Corrections were made based on the experimentally-estimated rheological model found with the multi-frequency torsional study.
Given the limited FOV in the experimental study, it was found that estimates of n or n 1 based on axially-or transversely-polarized (flexural) wave analysis were not reliable.In the finite element study, curve fits for these wavenumbers were based on wave motion 30 to 50 mm from the source, whereas in the experimental study the FOV only extended to 30 mm from the source.Signal-to-noise ratio (SNR) was too poor beyond that range.

V. DISCUSSION AND CONCLUDING REMARKS
The theoretical, numerical finite element and experimental studies of the previous sections have explored the confounding effects of finite dimensions and nonzero prestress on the elastography approach to estimating material viscoelastic properties in an isotropic cylindrical structure under uniaxial normal stress aligned with the cylinder axis.Additionally, a coordinate transformation approach-TAE-was introduced to estimate material viscoelastic properties independent of the prestress condition without requiring a priori knowledge of either the viscoelastic properties or stress conditions.Rather, only the amount of deformation, or strain, from the unstressed condition is required.Once viscoelastic properties are calculated, prestress can also be estimated.
The numerical and experimental studies show both the promise and implementation challenges of elastography and the TAE approach.In the numerical FEA study of Sec.III, the material shear storage l R and loss l I moduli, as well as the uniaxial normal stress r could be determined with accuracy within a few percent based on torsional-or axiallypolarized excitation, though accuracy generally degrades as the prestrain reaches 10% and 20%.The TAE approach formulated in Sec.II is inherently based on an assumption of linearity and as deformation increases accuracy will suffer.The approach could be improved by accounting for higher order terms in the quasi-static deformation model; this is left for future study.Related to this, in the present study it is assumed that the constitutive relation of elastic materials can be generalized to viscoelastic materials by substituting the real modulus with the complex modulus even if this linearization of the dynamic problem occurs after a finite static deformation (prestrain) has been applied.Further study is needed to evaluate this assumption.
A noted "divergence" between the FEA-based and experimental studies was the observation that the parameter l r was best approximated as À3/4 in the FEA simulation, whereas a value of -1/2 did better in the experimental study.Recall that this term accounts for the linear dependence of l on r.In the Refs.22 and 23 it is often denoted b 1 and is related to A used in other references, [24][25][26] where its value in soft tissue-like nearly incompressible materials can range between negative and positive values.It is hypothesized that its value depends on microstructure and thus can reveal material changes not captured by l. 24 Its approximate value of -1/2 in the experimental studies falls within the range found in the literature.The value of À3/4 found in the FEA studies, relative to -1/2, would lessen the effect of the prestress.This is commensurate with the fact that the linear material model used in the numerical study produces less of a change in stress levels for a given strain than does the Gent model, which has been shown to be accurate for the phantom material, Ecoflex-30. 33hereas the unprestressed torsional FEA studies perfectly reproduced the implemented rheological model over multiple frequencies, the multi-frequency experimental study illustrated the approximate nature of the rheological model under no prestress conditions.A least square error fit over the range from 200 to 600 Hz was used to estimate l 0 ; l a ; and a.This in turn was used in estimating the value of l r .If the baseline (unprestressed) rheological model is a poor fit, estimates of l r based on it would, of course, be inaccurate, as well.
Another source of error is the confounding effect of multiple wave types being present.This is particularly true for the axially-polarized and transverse-(flexural) polarized cases.Selecting a field of view away from the source and filtering can help with some of this, but not all of it.Even in the finite element studies, estimates of material shear properties based on flexural waves at any prestrain level were not accurate.
Simulations of a localized line segment source, an approximation of focused modulated radiation force of ultrasound that is commonly used in ultrasound-based elastography, highlighted how both small dimensions and prestress can alter estimates of the shear viscoelastic properties if those estimates are based on assuming bulk shear wave propagation and do not account for boundary and preloading effects.
Some "next steps" for advancing the strategy introduced here to decouple prestress and waveguide behavior from material shear stiffness estimates include consideration of more complex geometry and stress conditions, as well as anisotropic and nonuniform material properties.Finite element models based on medical images that can provide detailed geometry and localized deformation information under varying loading conditions and, if needed, measures of anisotropy and inhomogeneity, may provide a way to advance the TAE technique beyond simple geometries and assumptions of isotropy and homogeneity. 42

FIG. 1
FIG. 1. (Color online) Uniaxially prestressed cylinder with different harmonic excitation configurations.Shaded field of view regions FOV u ; FOV z ; and FOV x show the location of images presented in the following figures for torsionally-, axially-, and transversely-polarized wave motion.These are sagittal slices in the x-z or y-z plane.Deeper shaded blocks with vertical thick blue two-way arrows denote cross section of band surrounding the cylinder that delivers (z) axially-polarized harmonic actuation (b).Torsionally-polarized (a) and transversely-polarized (c) motion is input by appropriate motion at the base of the phantom, indicated by arrows.Axially-polarized (d) and transversely-polarized (e) line segment sources are also indicated at the geometric center of the phantom.
FIG. 5. (Color online) Shear loss versus shear storage modulus.The black dashed line is actual value based on Table I values for l 0 ; l a ; and a of the Fractional Voigt rheological model.Blue Â's are based on the torsionallypolarized curve fit to simulations at 500, 600, and 700 Hz with zero prestrain.Fitted values are also shown for prestrain of 5% without (red circle) and with (green asterisks) the TAE correction.

FIG. 6 .
FIG. 6. (Color online) Experimental setup for optical elastography using SLDV.Diagram with a bottom view and photo from the bottom and sideangle view showing cylindrical polymer phantom and tensioner.Colorcoded view from bottom: (a) Bore constraint cap; (b) piezo stack; (c) rotor hub; (d) image region of interest (ROI); (e) prestressed phantom; (f) tensioner inside phantom; (g) adjustable weight; (h) Polytec SLDV scan head.
FIG. 7. (Color online) Sagittal view of torsionally (u)-polarized in-phase wave motion from experimental study.(a)-(c) 0%, 10%, and 20% axial prestrain, respectively.Measurement differs from the FE result in Fig. 2 since the SLDV is measuring motion in the y-direction on the rounded surface of the cylindrical phantom, not at a central sagittal cross section.

TABLE I .
Parameter values for case studies.

TABLE II .
Wavenumber estimates for 35 mm phantom at 600 Hz at indicated prestrain levels.

TABLE III .
Wavenumber estimates for 20 mm phantom at 600 Hz at indicated prestrain levels.

TABLE IV .
Wavenumber estimates for 20 mm phantom at indicated prestrain levels.